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Abstract.  The  Green-Kubo  representation  of  the  shear  viscosity  of  a  dilute  gas  of  inelastic  hard  disks  derived  from  the 
Boltzmann  equation  is  evaluated  by  means  of  the  direct  simulation  Monte  Carlo  method.  The  relationship  between  the  one- 
particle  dynamics  of  the  original  expression  and  the  V-particle  dynamics  needed  for  the  simulation  is  analyzed.  The  presence 
of  velocity  correlations  in  the  latter  is  identified  and  their  possible  implications  are  discussed. 


INTRODUCTION 

The  prototypical  model  for  granular  gases  is  an  ensemble  of  smooth  inelastic  hard  spheres  or  disks  with  a  constant 
coefficient  of  normal  restitution  a.  For  this  model,  the  methods  of  kinetic  theory  have  been  applied  and,  in  particular, 
the  Boltzmann  equation  has  been  extended  [1,  2].  Then,  by  using  a  modification  of  the  Chapmann-Enskog  procedure, 
hydrodynamic  equations  have  been  derived  [3,  4].  As  the  reference  state,  the  isotropic  and  homogeneous  cooling  state 
(HCS)  of  a  granular  system  is  used.  Due  to  the  energy  dissipation,  the  HCS  is  not  stationary,  but  its  energy  decreases 
monotonically.  Hydrodynamics  has  been  employed  with  success  to  describe  the  behavior  of  granular  flows. 

As  for  molecular  systems,  the  Chapmann-Enskog  method  leads  to  complicated  differential  equations  that  have  to  be 
solved  in  some  approximations.  It  has  been  shown  that  the  aforementioned  equations  are  equivalent  to  a  representation 
of  the  transport  coefficients  that  has  a  similar  structure  to  the  Green-Kubo  relations  for  elastic  systems,  although  with 
significant  differences  [5].  Equivalent  expressions  have  been  derived  by  identifying  the  hydrodynamic  part  of  the 
spectrum  of  the  linearized  Boltzmann  equation  and  assuming  it  dominates  for  long  long  times  and  wavelengths  [6,  7]. 

Here,  an  evaluation  of  the  Green-Kubo  relation  for  the  shear  viscosity  by  means  of  A'-particlc  simulations,  namely 
the  direct  simulation  Monte  Carlo  (DSMC)  method  [8],  is  presented.  This  requires  going  from  the  one-particle 
effective  dynamics,  as  defined  by  a  linearized  Boltzmann  operator,  to  a  supposed  equivalent  A'-particlc  dynamics. 
The  relationship  between  both  descriptions  will  turn  out  to  be  important  for  the  analysis  of  the  results.  The  existence 
an  exact  mapping  of  the  HCS  onto  a  steady  state  will  be  exploited  as  the  basis  of  the  simulation  method  [9,  10]. 


THE  GREEN-KUBO  EXPRESSION  FOR  THE  SHEAR  VISCOSITY 


By  using  the  Chapman-Enskog  method  or  eigenfunctions  expansions  [5,  6],  the  shear  viscosity  r)  of  a  dilute  granular 
gas  of  N  smooth  inelastic  hard  disks  ( d  =  2)  or  spheres  (d  =  3)  of  mass  m  and  diameter  a  can  be  expressed  in  the  form 

t){T)  mnm£v0(T) J  dse~st=o/2  J  dci  XHCs(c\)A_xy(ci,s)^>2,xy(ci)-  (1) 

Here  n  is  the  density,  T  the  temperature,  vq(7')  =  (2k nT /in)1'2  with  kg  the  Boltzmann  constant,  i  =  (no'1  _1)_1 ,  and 


An  ( c )  —  C'xCy ,  d>2,.vy(c)  —  T 


dlnXncsjc) 

dcv 


(2) 


The  function  Xhcs{  Ci)  is  defined  from  the  solution  of  the  Boltzmann  equation  describing  the  HCS,  fncsis  i  •  f ),  as  [1] 


fHcsiy  id)  =  nv0d[T  (t)]xHCs{c\),  ci  = 


Vl 


voino]  ’ 


(3) 
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and  it  is  an  isotropic  function  of  the  vector  c  i .  Finally,  the  time  dependence  of  Axy  is  given  by 


A,v(c1^)=e'^(c^AAy(ci), 


(4) 


Ac(ci)  =  j dc2XHcs(c2)  J da ®{az  •  a)ci2-o[ba(ci,C2)  -  l](l  +  £gn)  + 


Co  d 
YCl'^7’ 


(5) 


where  C12  =  Ci  —  C2,  ©  is  the  Heaviside  step  function,  a  is  the  unit  vector  pointing  from  the  center  of  particle  2  to  the 
center  of  particle  1  at  contact,  the  operator  interchanges  the  subindexes  1  and  2  to  its  right,  and  ba(c\ . C2 )  is  an 
operator  replacing  all  the  velocities  c  1  and  C2  appearing  to  its  right  by  the  postcollisional  values  c'j  and  cC  given  by 


/  _ ,  1  +  a  ^  ^ 

C]  =  boC\  — Ci  —  (ct  •  Ci2jcr, 


/  ,  1  +  a.~ 

c2  see  ba c2  =  Cl  H - - — (a  ■  Ci2)a. 


(6) 


The  last  term  on  the  right  hand  side  of  Eq.  (5)  contains  the  time-independent  reduced  cooling  rate  of  the  HCS,  £0, 
debiting  the  evolution  of  the  temperature  in  that  state. 


d,T{t)  =  -i;HcS{T)T,  Co  = 


(’C HCS 


(7) 


vo(7T 

Equation  (1)  differs  from  the  standard  Green- Kubo  expression  for  dilute  molecular  systems  in  several  aspects: 

•  The  average  is  taken  over  the  velocity  distribution  of  the  HCS  and  not  over  the  Maxwellian  characterizing  the 
equilibrium  state.  Both  differ  for  a  <  1. 

•  The  time  correlation  function  couples  the  transversal  momentum  bux  A„  and  a  “modibed  flux”  <f>2 tXy. 

•  The  dynamics  includes  an  accelerating  streaming  between  collisions. 

•  The  time  integral  contains,  in  addition  to  the  correlation  function,  an  exponentially  decreasing  in  time  factor. 

To  render  Eq.  (1)  suitable  for  evaluation  by  means  of  /V-particle  simulations,  some  issues  must  be  addressed. 
First,  it  presents  the  technical  difficulty  that  Ac  involves  the  cooling  rate  Co  that,  therefore,  must  be  known  a  priori. 
Consequently,  it  is  useful  to  make  a  change  of  scale  [9,  10], 


_  Co 

1  =  ~ - s, 

2coo 


2coq£ 

W  =  —p — c 

Co 


with  C0q  being  an  arbitrary  time-independent  dimensionless  frequency.  Then,  Eq.  (1)  transforms  into 

nmv()(T ) 


*700  = 


vo,^7V 


J  die  a)°r (AAry(v,f)4>2rq>  {y/vo, st))st  ■ 


Here,  we  have  renamed  the  time  and  velocity  variables,  vo.,st  =  (2-knTst /in) 1  2  where 


fst  =  - -( 

and  the  angular  brackets  denote  average  debned  by 


"  “  kB  \  Co ;  ’ 


(a{y,t)b(\))st  =  dr  dy  fst(y)a(y,t)b(y),  fst{y)  =  nv0dXHCS  ( 

J  J  ’V  v’o,« 


(8) 


(9) 


(10) 


(ID 


Moreover,  the  time  dependence  is  now  given  by 


a(v,t)  =  erAs,('',')a{y ) 

with 

»  _  p  ^ 

A.sf(vi)  =  ad~1  J  dy2fs,(y2 )  J  da&(yl2  ■  a)yn  ■  a[ba(yuy2)  -  1](1  +  ^12)  +  fflovi  • 


(12) 

(13) 
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N-PARTICLE  REPRESENTATION 


In  Eq.  (9),  the  dynamics  is  defined  by  means  of  the  linear  operator  Ast.  Before  computating  it  by  /V-particle  simulation 
methods,  it  must  be  rewritten  in  terms  of  a  well  defined  particle  dynamics.  The  structure  of  Ast  suggests  to  introduce 
an  acceleration  streaming  between  collisions, 

^R/(r)  =  V/(f),  =*  fflDVi(f),  (14) 

while  the  effect  of  a  collision  between  particles  i  and  j  is  to  instantaneously  modify  their  velocities  accordingly  with 
the  same  rules  as  in  the  original  dynamics  and  given  in  Eqs.  (6),  i.e. 

v,  -  v'  =  baYi  =  Vi  -  ^  (a  ■  Yij)  a,  Yj  -  V'  ee  baYj  =  V,-  +  (5  •  Vy)  a.  (15) 

For  this  dynamics,  the  Liouville  equation  can  be  obtained.  Then,  by  means  of  any  of  the  standard  procedures,  the 
Boltzmann  equation  can  be  derived  in  the  low  density  limit.  It  has  the  form 

+  +  (16) 

where  /[/,/]  is  the  same  Boltzmann  collision  operator  as  in  the  original  dynamics.  It  is  now  easily  verified  that  fst  ( v ) 
is  a  steady  solution  of  this  equation  [10].  Let  us  assume  that  the  Liouville  equation  has  a  stationary  solution  pst(T), 
where  T  denotes  a  point  in  phase  space.  A  sufficient  condition  for  this  is  that  the  distribution  function  of  the  HCS  in 
the  actual,  original  dynamics  has  the  scaling  property 

pHCs(r,t)  =  (vo[Tm-dNPHCs({Kij,Vi/MTm)-  (17) 

Consider  time  correlation  functions  of  the  form 

Cab, st  (0  =  (A(t)B)  N,st  —  ( A)N,st(B)N,st ,  (18) 

where 

(A)N,st  =  J  dTpst(r)A(r),  (B)  N>st  =  J  dTpst(T)B(T), 

(A(t)B)N,sr  =  J  drps[(T)A(r,t)B(r), 

with 

N  N 

A(r):-£a(V,),  B(X)  =  '£b(\,). 

I- 1  1=1 

The  dynamical  variable  A  (T,f)  =  A  [L(r  )J  is  generated  from  A(T)  through  the  dynamics  defined  above.  The  low  density 
limit  of  Cab, st  can  be  obtained  by  using  the  same  procedures  as  for  molecular  systems  [10].  The  hypothesis  needed 
are  the  same  as  those  required  to  derive  the  Boltzmann  equation  [11],  plus  the  additional  assumption  that  velocity 
correlations  are  negligible  in  the  HCS  at  low  densities.  Then,  it  is  found  that  CAB,st  is  given  by  Eq.  (11).  This  provides 
a  direct  way  of  evaluating  the  latter  expression,  and  therefore  the  shear  viscosity,  by  means  of  /V-particlc  simulations, 
since  the  dynamics  defined  by  Eqs.  (14)  and  (15)  is  easily  implemented.  Of  course,  the  equality 

CAB,st  =  (a(\,t)b(\))s,  (22) 

only  holds  at  low  enough  densities  and  as  long  as  velocity  correlations  effects  are  negligible  in  the  HCS.  The  DSMC 
method  is  a  very  useful  tool  to  generate  the  low  density  dynamics  of  a  system  [8,  12].  This  method  is  designed  as  a 
real  /V-particle  dynamics  simulation  of  a  low  density  gas,  consequently  providing  its  complete  dynamical  description. 


(19) 

(20) 

(21) 
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FIGURE  1.  Plot  of  (lnX)'  =  <9  In Xhcs/^c  as  a  function  of  c  =  v/vo,s?  for  a  =  0.6.  The  circles  are  the  numerical  derivative  of 
the  simulation  results  and  the  solid  line  the  fitted  function.  The  velocity  is  measured  in  the  units  defined  in  the  text. 


SIMULATION  RESULTS 

We  have  performed  DSMC  simulations  of  a  system  of  N  =  104  inelastic  hard  disks.  Since  we  are  dealing  with  position 
independent  quantities  in  a  homogeneous  state,  it  is  enough  to  consider  one  cell  in  configuration  space.  This  allows  to 
increase  the  statistics  and  also  avoids  that  the  system  spontaneously  develops  spatial  inhomogeneities.  This  is  important 
since  the  HCS  is  unstable  with  respect  to  spatial  long  wavelength  perturbations  [13,  14], 

We  introduce  a  dimensionless  reduced  shear  viscosity  rj  *  by  scaling  the  viscosity  with  its  elastic  limit  in  the  first 
Sonine  approximation  rjq. 


V  = 


V(T) 

Vo(T)' 


Vo 


1  /  inksT  \ 

2cr  \  n  ) 


1/2 


Then,  Eq.  (9)  yields 


with 


V  = 


2\fm 

£v0,st 


[  dt  Jr,  (t)e~<0°t , 
Jo 


Jq  (t)  —  ^  (Ary  (  V .  f )<&2,xy  (  V  / Vo..,/  )  } st  • 

Let  us  remind  that  what  is  actually  computed  in  the  /V-particle  description  is 


(23) 


(24) 


(25) 


1  N  N 

Jq(j)  =  id)^2,xy(y j /v0,st)) N,st  ■ 


(26) 


The  simulations  show  that  the  system  reaches,  after  a  transient  period,  a  stationary  state  with  a  temperature  Tst.  Then, 
we  measured  the  velocity  function  d  In  Xhcs(c) / dc ,  c  =  v/vq.v,  to  determine  the  modified  flux  <t>2„ry(c),  defined  in  Eq. 
(2).  For  this,  the  range  of  v  was  partitioned  into  small  non-overlapping  bins  and  the  frequency  distribution  was  built. 
Afterwards,  the  numerical  derivative  was  computed.  An  example  of  the  results,  for  a  =  0.6,  is  given  in_Fig.  1.  Here 
and  in  the  following,  the  unit  of  mass  is  m,  the  unit  of  length  is  (  —  («(7)  ~ 1 ,  and  the  unit  of  time  is  (j2A:/i7’('0j//jzJ  1  2, 
where  T( 0)  is  the  initial  temperature.  Moreover,  we  set  kg  =  1,  implying  that  in  our  units  it  is  T( 0)  =  1/2.  For  large 
velocities,  d  In Xhcs(^)  j dc  tends  to  a  constant,  reflecting  the  exponential  decay  of  Xncs(c)  [15]. 

Once  the  function  <t>2..iy  has  been  determined  for  each  of  the  values  of  a  considered,  the  correlation  function  .//  (t  ) 
has  been  measured.  It  is  seen  to  decay  exponentially  in  time,  at  least  over  more  than  two  decades,  for  all  the  values  of 
the  restitution  coefficient  [16].  Then,  J'n  (t  )  was  fitted  to  an  exponential,  7/  (f)  =  7/  (O)e^  'LT,  and  the  fitted  parameters 
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FIGURE  2.  The  reduced  coefficient  of  shear  viscosity  r /  *  for  a  system  of  inelastic  disks  as  a  function  of  the  coefficient  of  normal 
restitution  a.  The  symbols  are  from  the  DSMC  method  and  the  solid  line  is  the  first  Sonine  approximation. 


JL  (0)  and  A,,  were  determined.  In  terms  of  them,  Eq.  (24)  reads 


2V27Zj'n(0) 
A,;  +  coq 


(27) 


The  values  of  r]*  obtained  in  this  way  are  shown  in  Fig.  2.  A  good  agreement  is  observed  with  the  results  obtained 
in  the  first  Sonine  approximation  [17],  although  a  slightly  different  behavior  is  also  clearly  identified.  It  is  interesting 
to  consider  the  initial  value  of  JL  (t  )  and  decompose  it  in  the  form 


j;(0)=y;(1)(0)+42)(0), 


(28) 


where 


’(0) 


1  N 

IL(Avv(V;)‘I>2^(v//vo,sf))iV,v  , 


42)(o) 


1  N  N 

tj  L  L(A^(v0^2 sy(yj/vo,st))N,st  ■ 
i  Mi 


(29) 


The  first  component,  fl 1 1  (0)  is  identically  the  same  as  Jn  (0  ),  where  Jn  (t)  was  defined  in  Eq.  (25).  It  is  independent 
of  the  velocity  correlations  present  in  the  system.  In  fact,  its  exact  value  can  be  easily  obtained,  and  for  the  choice 
of  o\)  used  in  the  simulations  and  the  units  we  are  employing,  it  is  4*^(0)  =  Jt\ ( 0 )  =  1/2  for  all  a  [16].  The  other 

(2) 

part,  Jn  ( 0),  contains  the  effect  of  velocity  correlations  and  was  assumed  to  be  negligible.  The  simulation  results  for 

JL  (0)  are  plotted  in  Fig.  3.  For  a  <  0.6,  contributions  from  J,j  become  significant,  indicating  the  presence  of  velocity 
correlations  in  the  system.  This  result  raises  some  fundamental  questions.  Are  the  observed  correlations  an  artifact  of 
the  DSMC  method?  We  believe  they  are  not,  since  velocity  correlations  in  the  HCS  also  have  been  recently  observed 
and  measured  in  molecular  dynamics  simulations  [18].  Why  there  is  such  a  good  agreement  between  the  first  Sonine 
approximation  predictions  and  the  simulation  results  for  the  shear  viscosity,  even  from  strong  dissipation?  A  possible 
answer  is  that  the  Sonine  expansion  is  introduced,  and  truncated,  both  in  the  initial  conditions  and  in  the  dynamics  of 
the  fluxes,  leading  to  some  kind  of  compensation.  For  small  values  of  a,  which  Green-Kubo  expression  is  the  right 
one,  that  with  Jr\{t)  or  that  with  .7/  (f)?,  or  are  both  wrong?  And,  finally,  can  it  be  concluded  from  these  results  that 
the  Boltzmann  equation  should  be  somehow  revised  for  strong  inelasticity?  All  these  points  clearly  deserve  further 
analysis. 

Although  in  this  paper  we  have  restricted  ourselves  to  the  shear  viscosity,  the  Green-Kubo  expressions  for  the 
other  two  transport  coefficients,  associated  with  the  heat  flux,  can  be  computed  in  a  similar  way.  For  not  too  strong 
dissipation,  the  simulation  results  agree  quite  well  with  the  first  Sonine  approximation  predictions.  Nevertheless,  the 
discrepancies  grow  rather  fast  as  a  decreases  below  a  ~  0.7.  This  might  be,  at  least  partially,  due  to  the  presence 
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FIGURE  3.  DSMC  results  for  the  correlation  function  J'^  (0)  (filled  circles)  and  its  diagonal  part  0)  (empty  circles). 


of  velocity  correlations  in  the  DSMC  results.  A  detailed  discussion  of  these  transport  coefficients  will  be  published 
elsewhere  [16]. 
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